#---- Replication script 3/4  ----

#Turnbull-Dugarte
#"Liberal in the sheets [...]"

###This script produces all analysis from Study II

remove(list=ls())

#install.packages("pacman")
library(pacman)

p_load(here, tidyverse, dplyr, haven, modelsummary, rcartcolor, jtools, broom, sdamr,
       distributional, list, ggpubr, sruvey, xtable, interactions, viridisLite,
       ggdist, patchwork, estimatr, margins, grf, starbility, miceadds, modelr)

colors <- viridis(option = "plasma", 
                  begin = 0.2, 
                  end = 0.8, 
                  direction = -1, 
                  n = 4)


##Loading data##
df_obs <- read_dta("data/study1_data_observational.dta")
df_obs <-as.data.frame(df_obs)

df_het_obs <- df_obs %>% 
  filter(sexuality==0)
df_het_obs <- df_het_obs %>% 
  filter(gender!=2)

df_het_obs<- df_het_obs  %>% 
  mutate(countryID=as.factor(countryID),
         gender=as.factor(gender),
         degree=as.factor(degree),
         age=as.numeric(age))

cbp2 <- c( "#FC4C02", "purple")

mod_full<- lm(direct_bi ~ gender, data = df_het_obs)
df_het_obs$predicted<-predict(mod_full, df_het_obs)
summary(mod_full)
table(df_het_obs$predicted)


plot1<- df_het_obs %>%
  data_grid(gender, .model = mod_full) %>%
  augment(mod_full, newdata = ., se_fit = TRUE) %>%
  ggplot(aes(x = gender, fill=gender)) +
  scale_fill_manual(values = cbp2) +
  scale_color_manual(values = cbp2)+
  geom_jitter(data=df_het_obs, aes(x=gender, y=predicted, fill=gender),
              height=.25, width=.2,size=3, alpha=.1,
              pch=21, color="grey57")+  
  stat_dist_halfeye(
    aes(dist = dist_student_t(df = df.residual(mod_full), mu = .fitted, sigma = .se.fit)), 
    scale = .3, position = position_nudge(x = .35))+
  geom_hline(yintercept=.5, linetype="dashed")+
  labs(y="", x="",
       caption = "Confidence intervals at 95% (N:3087)",
       title="Self-reported prejudice against bisexual partners",
       subtitle="Pr(Considers being bisexual a 'red flag')")+
  guides(color=FALSE)+
  ylim(0,1)+
  theme_ggdist()+
  theme(legend.position = "none",
        plot.title = element_text(color="black", face=2, hjust=0),
        axis.text.x = element_text(face=2),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank())+
  coord_flip()+
  annotate(
    geom="label", x = 1, y = 0.42, size = 3, colour="#FC4C02", fontface=2,
    label = "0.42", hjust=0.5)+
  annotate(
    geom="label", x = 2, y = 0.58, size = 3, colour="purple", fontface=2,
    label = "0.58", hjust=0.5)+
  annotate(
    geom="label", x = 1, y = .1, size = 3, colour="#FC4C02", fontface=2,
    label = "Men", hjust=0.5)+
  annotate(
    geom="label", x = 2, y = .1, size = 3, colour="purple", fontface=2,
    label = "Women", hjust=0.5)+
  annotate("segment", x = 1, xend = 2, y = .9, yend = .9, color = "black", linetype = "solid")+
  annotate("segment", x = 2, xend = 2, y = .88, yend = .9, color = "black", linetype = "solid")+
  annotate("segment", x = 1, xend = 1, y = .88, yend = .9, color = "black", linetype = "solid")+
  annotate("label", x = 1.5, y = .9, label = "p<0.001", size = 2, color = "black")
plot1
ggsave(filename = here("figures/Figure3.png"),
       height=4, dpi = 600)
ggsave(filename = here("figures/tifs/Figure3.tif"),
       height=4, dpi = 600)



##Loading data for List##
df_list <- read_dta("data/study1_data_list.dta")
df_list <-as.data.frame(df_list)
df_list<- df_list  %>% 
  mutate(treat=as.factor(treat),
         IDvar=as.factor(IDvar),
         countryID=as.factor(countryID),
         gender=as.factor(gender),
         round=as.factor(round),
         list=as.factor(list))

df_het <- df_list %>% 
  filter(sexuality==0)
df_men <- df_het %>% 
  filter(gender==0)
df_women <- df_het %>% 
  filter(gender==1)
df_hetUK <- df_het %>% 
  filter(country=="UK")
df_menUK <- df_hetUK  %>% 
  filter(gender==0)
df_womenUK <- df_hetUK  %>% 
  filter(gender==1)
df_hetES <- df_het %>% 
  filter(country=="ES")
df_menES <- df_hetES  %>% 
  filter(gender==0)
df_womenES <- df_hetES  %>% 
  filter(gender==1)



mod_full<- lm_robust(list_answer ~ treat + round + list, data = df_het, clusters = IDvar)
tidy_full <- tidy(mod_full)
tidy_full$sample <- "Full sample"
tidy_full$outcome <- "Peers"

mod_full2<- lm_robust(list_answer_self ~ treat + round + list, data = df_het, clusters = IDvar)
tidy_full2 <- tidy(mod_full2)
tidy_full2$sample <- "Full sample"
tidy_full2$outcome <- "Self"

mod_men<- lm_robust(list_answer ~ treat + round + list, data = df_men, clusters = IDvar)
tidy_men <- tidy(mod_men)
tidy_men$sample <- "Men"
tidy_men$outcome <- "Peers"

mod_men2<- lm_robust(list_answer_self ~ treat + round + list, data = df_men, clusters = IDvar)
tidy_men2 <- tidy(mod_men2)
tidy_men2$sample <- "Men"
tidy_men2$outcome <- "Self"

mod_women<- lm_robust(list_answer ~ treat + round + list, data = df_women, clusters = IDvar)
tidy_women <- tidy(mod_women)
tidy_women$sample <- "Women"
tidy_women$outcome <- "Peers"

mod_women2<- lm_robust(list_answer_self ~ treat + round + list, data = df_women, clusters = IDvar)
tidy_women2 <- tidy(mod_women2)
tidy_women2$sample <- "Women"
tidy_women2$outcome <- "Self"

models<- rbind(tidy_full, tidy_men, tidy_women, tidy_full2, tidy_men2, tidy_women2)
models<-subset(models, term %in% c("treat1")) 

models$lower90 <- models$estimate - (1.645 * models$std.error)
models$higher90 <- models$estimate + (1.645 * models$std.error)
models$lower95 <- models$estimate - (1.96 * models$std.error)
models$higher95 <- models$estimate + (1.96 * models$std.error)
models$lower99 <- models$estimate - (2.58 * models$std.error)
models$higher99 <- models$estimate + (2.58 * models$std.error)

col2<- c("steelblue", "#B12A90FF")

plot_full<- ggplot(models, aes(y = sample, x=estimate, color=outcome)) + 
  geom_vline(xintercept = 0, linetype="longdash", colour="black", alpha=.5)+
  scale_color_manual(values = col2, guide = "none")+
  xlim(0,.8)+
  geom_segment(data = subset(models, outcome == "Peers"),
               aes(x = lower90, xend = higher90, y = sample, yend = sample),
               size = 4, alpha = 1,  position = position_nudge(y = .1)) + 
  geom_segment(data = subset(models, outcome == "Peers"),
               aes(x = lower95, xend = higher95, y = sample, yend = sample),
               size = 4, alpha = .6,  position = position_nudge(y = .1)) + 
  geom_segment(data = subset(models, outcome == "Peers"),
               aes(x = lower99, xend = higher99, y = sample, yend = sample),
               size = 4, alpha = .3, position = position_nudge(y = .1)) + 
  geom_point(data = subset(models, outcome == "Peers"),
             aes(x = estimate, y = sample), fill = "white", color = "grey67", size = 2.5, pch = 21,  
             position = position_nudge(y = .1))+
  geom_segment(data = subset(models, outcome == "Self"),
               aes(x = lower90, xend = higher90, y = sample, yend = sample),
               size = 4, alpha = 1,  position = position_nudge(y = -.1)) + 
  geom_segment(data = subset(models, outcome == "Self"),
               aes(x = lower95, xend = higher95, y = sample, yend = sample),
               size = 4, alpha = .6,  position = position_nudge(y = -.1)) + 
  geom_segment(data = subset(models, outcome == "Self"),
               aes(x = lower99, xend = higher99, y = sample, yend = sample),
               size = 4, alpha = .3, position = position_nudge(y = -.1)) + 
  geom_point(data = subset(models, outcome == "Self"),
             aes(x = estimate, y = sample), fill = "white", color = "grey67", size = 2.5, pch = 21,  
             position = position_nudge(y = -.1))+
  geom_text(data = subset(models, outcome == "Self"),aes(label=sprintf("%.2f", 
             round(estimate, digits = 2))), nudge_y = -.24, size=2.5)+
  geom_text(data = subset(models, outcome == "Peers"),aes(label=sprintf("%.2f", 
             round(estimate, digits = 2))), nudge_y = .24, size=2.5)+
  theme_ggdist()+
  labs(x="",
       title="Combined sample")+
  theme(legend.position = "none",
        axis.title.y = element_blank(),
        plot.title=element_text(color="black", face="bold", hjust=.5),
        axis.text.x=element_text(color="black"),
        axis.text.y = element_text(color="black", face="bold"))


mod_fullUK<- lm_robust(list_answer ~ treat + round + list, data = df_hetUK)
tidy_fullUK <- tidy(mod_fullUK)
tidy_fullUK$sample <- "Full sample"
tidy_fullUK$outcome <- "Peers"

mod_fullUK2<- lm_robust(list_answer_self ~ treat + round + list, data = df_hetUK)
tidy_fullUK2 <- tidy(mod_fullUK2)
tidy_fullUK2$sample <- "Full sample"
tidy_fullUK2$outcome <- "Self"

mod_menUK<- lm_robust(list_answer ~ treat + round + list, data = df_menUK)
tidy_menUK <- tidy(mod_menUK)
tidy_menUK$sample <- "Men"
tidy_menUK$outcome <- "Peers"

mod_menUK2<- lm_robust(list_answer_self ~ treat + round + list, data = df_menUK)
tidy_menUK2 <- tidy(mod_menUK2)
tidy_menUK2$sample <- "Men"
tidy_menUK2$outcome <- "Self"

mod_womenUK<- lm_robust(list_answer ~ treat + round + list, data = df_womenUK)
tidy_womenUK <- tidy(mod_womenUK)
tidy_womenUK$sample <- "Women"
tidy_womenUK$outcome <- "Peers"

mod_womenUK2<- lm_robust(list_answer_self ~ treat + round + list, data = df_womenUK)
tidy_womenUK2 <- tidy(mod_womenUK2)
tidy_womenUK2$sample <- "Women"
tidy_womenUK2$outcome <- "Self"

modelsUK<- rbind(tidy_fullUK, tidy_menUK, tidy_womenUK, tidy_fullUK2, tidy_menUK2, tidy_womenUK2)
modelsUK<-subset(modelsUK, term %in% c("treat1")) 

modelsUK$lower90 <- modelsUK$estimate - (1.645 * modelsUK$std.error)
modelsUK$higher90 <- modelsUK$estimate + (1.645 * modelsUK$std.error)
modelsUK$lower95 <- modelsUK$estimate - (1.96 * modelsUK$std.error)
modelsUK$higher95 <- modelsUK$estimate + (1.96 * modelsUK$std.error)
modelsUK$lower99 <- modelsUK$estimate - (2.58 * modelsUK$std.error)
modelsUK$higher99 <- modelsUK$estimate + (2.58 * modelsUK$std.error)



plot_UK<- ggplot(modelsUK, aes(y = sample, x=estimate, color=outcome)) + 
  geom_vline(xintercept = 0, linetype="longdash", colour="black", alpha=.5)+
  scale_color_manual(values = col2, guide = "none")+
  xlim(0,.8)+
  geom_segment(data = subset(modelsUK, outcome == "Peers"),
               aes(x = lower90, xend = higher90, y = sample, yend = sample),
               size = 4, alpha = 1,  position = position_nudge(y = .1)) + 
  geom_segment(data = subset(modelsUK, outcome == "Peers"),
               aes(x = lower95, xend = higher95, y = sample, yend = sample),
               size = 4, alpha = .6,  position = position_nudge(y = .1)) + 
  geom_segment(data = subset(modelsUK, outcome == "Peers"),
               aes(x = lower99, xend = higher99, y = sample, yend = sample),
               size = 4, alpha = .3, position = position_nudge(y = .1)) + 
  geom_point(data = subset(modelsUK, outcome == "Peers"),
             aes(x = estimate, y = sample), fill = "white", color = "grey67", size = 2.5, pch = 21,  
             position = position_nudge(y = .1))+
  geom_segment(data = subset(modelsUK, outcome == "Self"),
               aes(x = lower90, xend = higher90, y = sample, yend = sample),
               size = 4, alpha = 1,  position = position_nudge(y = -.1)) + 
  geom_segment(data = subset(modelsUK, outcome == "Self"),
               aes(x = lower95, xend = higher95, y = sample, yend = sample),
               size = 4, alpha = .6,  position = position_nudge(y = -.1)) + 
  geom_segment(data = subset(modelsUK, outcome == "Self"),
               aes(x = lower99, xend = higher99, y = sample, yend = sample),
               size = 4, alpha = .3, position = position_nudge(y = -.1)) + 
  geom_point(data = subset(modelsUK, outcome == "Self"),
             aes(x = estimate, y = sample), fill = "white", color = "grey67", size = 2.5, pch = 21,  
             position = position_nudge(y = -.1))+
  geom_text(data = subset(modelsUK, outcome == "Self"),aes(label=sprintf("%.2f", 
                                                                          round(estimate, digits = 2))), nudge_y = -.24, size=2.5)+
  geom_text(data = subset(modelsUK, outcome == "Peers"),aes(label=sprintf("%.2f", 
                                                                           round(estimate, digits = 2))), nudge_y = .24, size=2.5)+
  theme_ggdist()+
  annotate("label", x = 0, y = 3.4, label = "Red flag for peers",
           color="steelblue", size=3, fontface="bold", hjust=0)+
  annotate("label", x = 0, y = 3.2, label = "Red flag for self",
           color="#B12A90FF", size=3, fontface="bold", hjust=0)+
  labs(x="",
       title="Britain")+
  theme(legend.position = "none",
        axis.title.y = element_blank(),
        plot.title=element_text(color="black", face="bold", hjust=.5),
        axis.text.x=element_text(color="black"),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank())


mod_fullES<- lm_robust(list_answer ~ treat + round + list, data = df_hetES)
tidy_fullES <- tidy(mod_fullES)
tidy_fullES$sample <- "Full sample"
tidy_fullES$outcome <- "Peers"

mod_fullES2<- lm_robust(list_answer_self ~ treat + round + list, data = df_hetES)
tidy_fullES2 <- tidy(mod_fullES2)
tidy_fullES2$sample <- "Full sample"
tidy_fullES2$outcome <- "Self"

mod_menES<- lm_robust(list_answer ~ treat + round + list, data = df_menES)
tidy_menES <- tidy(mod_menES)
tidy_menES$sample <- "Men"
tidy_menES$outcome <- "Peers"

mod_menES2<- lm_robust(list_answer_self ~ treat + round + list, data = df_menES)
tidy_menES2 <- tidy(mod_menES2)
tidy_menES2$sample <- "Men"
tidy_menES2$outcome <- "Self"

mod_womenES<- lm_robust(list_answer ~ treat + round + list, data = df_womenES)
tidy_womenES <- tidy(mod_womenES)
tidy_womenES$sample <- "Women"
tidy_womenES$outcome <- "Peers"

mod_womenES2<- lm_robust(list_answer_self ~ treat + round + list, data = df_womenES)
tidy_womenES2 <- tidy(mod_womenES2)
tidy_womenES2$sample <- "Women"
tidy_womenES2$outcome <- "Self"

modelsES<- rbind(tidy_fullES, tidy_menES, tidy_womenES, tidy_fullES2, tidy_menES2, tidy_womenES2)
modelsES<-subset(modelsES, term %in% c("treat1")) 

modelsES$lower90 <- modelsES$estimate - (1.645 * modelsES$std.error)
modelsES$higher90 <- modelsES$estimate + (1.645 * modelsES$std.error)
modelsES$lower95 <- modelsES$estimate - (1.96 * modelsES$std.error)
modelsES$higher95 <- modelsES$estimate + (1.96 * modelsES$std.error)
modelsES$lower99 <- modelsES$estimate - (2.58 * modelsES$std.error)
modelsES$higher99 <- modelsES$estimate + (2.58 * modelsES$std.error)



plot_ES<- ggplot(modelsES, aes(y = sample, x=estimate, color=outcome)) + 
  geom_vline(xintercept = 0, linetype="longdash", colour="black", alpha=.5)+
  scale_color_manual(values = col2, guide = "none")+
  xlim(0,.8)+
  geom_segment(data = subset(modelsES, outcome == "Peers"),
               aes(x = lower90, xend = higher90, y = sample, yend = sample),
               size = 4, alpha = 1,  position = position_nudge(y = .1)) + 
  geom_segment(data = subset(modelsES, outcome == "Peers"),
               aes(x = lower95, xend = higher95, y = sample, yend = sample),
               size = 4, alpha = .6,  position = position_nudge(y = .1)) + 
  geom_segment(data = subset(modelsES, outcome == "Peers"),
               aes(x = lower99, xend = higher99, y = sample, yend = sample),
               size = 4, alpha = .3, position = position_nudge(y = .1)) + 
  geom_point(data = subset(modelsES, outcome == "Peers"),
             aes(x = estimate, y = sample), fill = "white", color = "grey67", size = 2.5, pch = 21,  
             position = position_nudge(y = .1))+
  geom_segment(data = subset(modelsES, outcome == "Self"),
               aes(x = lower90, xend = higher90, y = sample, yend = sample),
               size = 4, alpha = 1,  position = position_nudge(y = -.1)) + 
  geom_segment(data = subset(modelsES, outcome == "Self"),
               aes(x = lower95, xend = higher95, y = sample, yend = sample),
               size = 4, alpha = .6,  position = position_nudge(y = -.1)) + 
  geom_segment(data = subset(modelsES, outcome == "Self"),
               aes(x = lower99, xend = higher99, y = sample, yend = sample),
               size = 4, alpha = .3, position = position_nudge(y = -.1)) + 
  geom_point(data = subset(modelsES, outcome == "Self"),
             aes(x = estimate, y = sample), fill = "white", color = "grey67", size = 2.5, pch = 21,  
             position = position_nudge(y = -.1))+
  geom_text(data = subset(modelsES, outcome == "Self"),aes(label=sprintf("%.2f", 
                                                                          round(estimate, digits = 2))), nudge_y = -.24, size=2.5)+
  geom_text(data = subset(modelsES, outcome == "Peers"),aes(label=sprintf("%.2f", 
                                                                           round(estimate, digits = 2))), nudge_y = .24, size=2.5)+
  theme_ggdist()+
  labs(x="",
       title="Spain")+
  theme(legend.position = "none",
        axis.title.y = element_blank(),
        plot.title=element_text(color="black", face="bold", hjust=.5),
        axis.text.x=element_text(color="black"),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank())


mod_diff_full<- lm_robust(list_answer ~ treat*gender + round + list, data = df_het, clusters = IDvar)
tidy_diff_full <- tidy(mod_diff_full)
tidy_diff_full$sample <- "Full sample"
tidy_diff_full$outcome <- "Peers"

mod_diff_full2<- lm_robust(list_answer_self ~ treat*gender + round + list, data = df_het, clusters = IDvar)
tidy_diff_full2 <- tidy(mod_diff_full2)
tidy_diff_full2$sample <- "Full sample"
tidy_diff_full2$outcome <- "Self"

models_diff<- rbind(tidy_diff_full, tidy_diff_full2)
models_diff<-subset(models_diff, term %in% c("treat1:gender1")) 
models_diff <- models_diff %>%
  mutate(sample= ifelse(sample == "Full sample", "Gender diff.\nin ATE", sample))

models_diff$lower90 <- models_diff$estimate - (1.645 * models_diff$std.error)
models_diff$higher90 <- models_diff$estimate + (1.645 * models_diff$std.error)
models_diff$lower95 <- models_diff$estimate - (1.96 * models_diff$std.error)
models_diff$higher95 <- models_diff$estimate + (1.96 * models_diff$std.error)
models_diff$lower99 <- models_diff$estimate - (2.58 * models_diff$std.error)
models_diff$higher99 <- models_diff$estimate + (2.58 * models_diff$std.error)

col2<- c("steelblue", "#B12A90FF")

diff_full_plot<- ggplot(models_diff, aes(y = sample, x=estimate, color=outcome)) + 
  geom_vline(xintercept = 0, linetype="longdash", colour="black", alpha=.5)+
  scale_color_manual(values = col2, guide = "none")+
  xlim(-.25,.45)+
  geom_segment(data = subset(models_diff, outcome == "Peers"),
               aes(x = lower90, xend = higher90, y = sample, yend = sample),
               size = 4, alpha = 1,  position = position_nudge(y = .2)) + 
  geom_segment(data = subset(models_diff, outcome == "Peers"),
               aes(x = lower95, xend = higher95, y = sample, yend = sample),
               size = 4, alpha = .6,  position = position_nudge(y = .2)) + 
  geom_segment(data = subset(models_diff, outcome == "Peers"),
               aes(x = lower99, xend = higher99, y = sample, yend = sample),
               size = 4, alpha = .3, position = position_nudge(y = .2)) + 
  geom_point(data = subset(models_diff, outcome == "Peers"),
             aes(x = estimate, y = sample), fill = "white", color = "grey67", size = 2.5, pch = 21,  
             position = position_nudge(y = .2))+
  geom_segment(data = subset(models_diff, outcome == "Self"),
               aes(x = lower90, xend = higher90, y = sample, yend = sample),
               size = 4, alpha = 1,  position = position_nudge(y = -.2)) + 
  geom_segment(data = subset(models_diff, outcome == "Self"),
               aes(x = lower95, xend = higher95, y = sample, yend = sample),
               size = 4, alpha = .6,  position = position_nudge(y = -.2)) + 
  geom_segment(data = subset(models_diff, outcome == "Self"),
               aes(x = lower99, xend = higher99, y = sample, yend = sample),
               size = 4, alpha = .3, position = position_nudge(y = -.2)) + 
  geom_point(data = subset(models_diff, outcome == "Self"),
             aes(x = estimate, y = sample), fill = "white", color = "grey67", size = 2.5, pch = 21,  
             position = position_nudge(y = -.2))+
  geom_text(data = subset(models_diff, outcome == "Self"),aes(label=sprintf("%.2f", 
                                                                            round(estimate, digits = 2))), nudge_y = -.35, size=2.5)+
  geom_text(data = subset(models_diff, outcome == "Peers"),aes(label=sprintf("%.2f", 
                                                                             round(estimate, digits = 2))), nudge_y = .35, size=2.5)+
  theme_ggdist()+
  labs(x="",
       title="")+
  theme(legend.position = "none",
        axis.title.y = element_blank(),
        plot.title=element_text(color="black", face="bold", hjust=.5),
        axis.text.x=element_text(color="black"),
        axis.text.y = element_text(color="black", face="bold"))

mod_diff_ES<- lm_robust(list_answer ~ treat*gender + round + list, data = df_hetES, clusters = IDvar)
tidy_diff_ES <- tidy(mod_diff_ES)
tidy_diff_ES$sample <- "Full sample"
tidy_diff_ES$outcome <- "Peers"

mod_diff_ES2<- lm_robust(list_answer_self ~ treat*gender + round + list, data = df_hetES, clusters = IDvar)
tidy_diff_ES2 <- tidy(mod_diff_ES2)
tidy_diff_ES2$sample <- "Full sample"
tidy_diff_ES2$outcome <- "Self"

models_diff_ES<- rbind(tidy_diff_ES, tidy_diff_ES2)
models_diff_ES<-subset(models_diff_ES, term %in% c("treat1:gender1")) 

models_diff_ES$lower90 <- models_diff_ES$estimate - (1.645 * models_diff_ES$std.error)
models_diff_ES$higher90 <- models_diff_ES$estimate + (1.645 * models_diff_ES$std.error)
models_diff_ES$lower95 <- models_diff_ES$estimate - (1.96 * models_diff_ES$std.error)
models_diff_ES$higher95 <- models_diff_ES$estimate + (1.96 * models_diff_ES$std.error)
models_diff_ES$lower99 <- models_diff_ES$estimate - (2.58 * models_diff_ES$std.error)
models_diff_ES$higher99 <- models_diff_ES$estimate + (2.58 * models_diff_ES$std.error)

col2<- c("steelblue", "#B12A90FF")

diffES_plot<- ggplot(models_diff_ES, aes(y = sample, x=estimate, color=outcome)) + 
  geom_vline(xintercept = 0, linetype="longdash", colour="black", alpha=.5)+
  scale_color_manual(values = col2, guide = "none")+
  xlim(-.25,.45)+
  geom_segment(data = subset(models_diff_ES, outcome == "Peers"),
               aes(x = lower90, xend = higher90, y = sample, yend = sample),
               size = 4, alpha = 1,  position = position_nudge(y = .2)) + 
  geom_segment(data = subset(models_diff_ES, outcome == "Peers"),
               aes(x = lower95, xend = higher95, y = sample, yend = sample),
               size = 4, alpha = .6,  position = position_nudge(y = .2)) + 
  geom_segment(data = subset(models_diff_ES, outcome == "Peers"),
               aes(x = lower99, xend = higher99, y = sample, yend = sample),
               size = 4, alpha = .3, position = position_nudge(y = .2)) + 
  geom_point(data = subset(models_diff_ES, outcome == "Peers"),
             aes(x = estimate, y = sample), fill = "white", color = "grey67", size = 2.5, pch = 21,  
             position = position_nudge(y = .2))+
  geom_segment(data = subset(models_diff_ES, outcome == "Self"),
               aes(x = lower90, xend = higher90, y = sample, yend = sample),
               size = 4, alpha = 1,  position = position_nudge(y = -.2)) + 
  geom_segment(data = subset(models_diff_ES, outcome == "Self"),
               aes(x = lower95, xend = higher95, y = sample, yend = sample),
               size = 4, alpha = .6,  position = position_nudge(y = -.2)) + 
  geom_segment(data = subset(models_diff_ES, outcome == "Self"),
               aes(x = lower99, xend = higher99, y = sample, yend = sample),
               size = 4, alpha = .3, position = position_nudge(y = -.2)) + 
  geom_point(data = subset(models_diff_ES, outcome == "Self"),
             aes(x = estimate, y = sample), fill = "white", color = "grey67", size = 2.5, pch = 21,  
             position = position_nudge(y = -.2))+
  geom_text(data = subset(models_diff_ES, outcome == "Self"),aes(label=sprintf("%.2f", 
                                                                               round(estimate, digits = 2))), nudge_y = -.35, size=2.5)+
  geom_text(data = subset(models_diff_ES, outcome == "Peers"),aes(label=sprintf("%.2f", 
                                                                                round(estimate, digits = 2))), nudge_y = .35, size=2.5)+
  theme_ggdist()+
  labs(x="",
       title="")+
  theme(legend.position = "none",
        axis.title.y = element_blank(),
        plot.title=element_text(color="black", face="bold", hjust=.5),
        axis.text.x=element_text(color="black"),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank())


mod_diff_UK<- lm_robust(list_answer ~ treat*gender + round + list, data = df_hetUK, clusters = IDvar)
tidy_diff_UK <- tidy(mod_diff_UK)
tidy_diff_UK$sample <- "Full sample"
tidy_diff_UK$outcome <- "Peers"

mod_diff_UK2<- lm_robust(list_answer_self ~ treat*gender + round + list, data = df_hetUK, clusters = IDvar)
tidy_diff_UK2 <- tidy(mod_diff_UK2)
tidy_diff_UK2$sample <- "Full sample"
tidy_diff_UK2$outcome <- "Self"

models_diff_UK<- rbind(tidy_diff_UK, tidy_diff_UK2)
models_diff_UK<-subset(models_diff_UK, term %in% c("treat1:gender1")) 

models_diff_UK$lower90 <- models_diff_UK$estimate - (1.645 * models_diff_UK$std.error)
models_diff_UK$higher90 <- models_diff_UK$estimate + (1.645 * models_diff_UK$std.error)
models_diff_UK$lower95 <- models_diff_UK$estimate - (1.96 * models_diff_UK$std.error)
models_diff_UK$higher95 <- models_diff_UK$estimate + (1.96 * models_diff_UK$std.error)
models_diff_UK$lower99 <- models_diff_UK$estimate - (2.58 * models_diff_UK$std.error)
models_diff_UK$higher99 <- models_diff_UK$estimate + (2.58 * models_diff_UK$std.error)

col2<- c("steelblue", "#B12A90FF")

diffUK_plot<- ggplot(models_diff_UK, aes(y = sample, x=estimate, color=outcome)) + 
  geom_vline(xintercept = 0, linetype="longdash", colour="black", alpha=.5)+
  scale_color_manual(values = col2, guide = "none")+
  xlim(-.25,.45)+
  geom_segment(data = subset(models_diff_UK, outcome == "Peers"),
               aes(x = lower90, xend = higher90, y = sample, yend = sample),
               size = 4, alpha = 1,  position = position_nudge(y = .2)) + 
  geom_segment(data = subset(models_diff_UK, outcome == "Peers"),
               aes(x = lower95, xend = higher95, y = sample, yend = sample),
               size = 4, alpha = .6,  position = position_nudge(y = .2)) + 
  geom_segment(data = subset(models_diff_UK, outcome == "Peers"),
               aes(x = lower99, xend = higher99, y = sample, yend = sample),
               size = 4, alpha = .3, position = position_nudge(y = .2)) + 
  geom_point(data = subset(models_diff_UK, outcome == "Peers"),
             aes(x = estimate, y = sample), fill = "white", color = "grey67", size = 2.5, pch = 21,  
             position = position_nudge(y = .2))+
  geom_segment(data = subset(models_diff_UK, outcome == "Self"),
               aes(x = lower90, xend = higher90, y = sample, yend = sample),
               size = 4, alpha = 1,  position = position_nudge(y = -.2)) + 
  geom_segment(data = subset(models_diff_UK, outcome == "Self"),
               aes(x = lower95, xend = higher95, y = sample, yend = sample),
               size = 4, alpha = .6,  position = position_nudge(y = -.2)) + 
  geom_segment(data = subset(models_diff_UK, outcome == "Self"),
               aes(x = lower99, xend = higher99, y = sample, yend = sample),
               size = 4, alpha = .3, position = position_nudge(y = -.2)) + 
  geom_point(data = subset(models_diff_UK, outcome == "Self"),
             aes(x = estimate, y = sample), fill = "white", color = "grey67", size = 2.5, pch = 21,  
             position = position_nudge(y = -.2))+
  geom_text(data = subset(models_diff_UK, outcome == "Self"),aes(label=sprintf("%.2f", 
                                                                               round(estimate, digits = 2))), nudge_y = -.35, size=2.5)+
  geom_text(data = subset(models_diff_UK, outcome == "Peers"),aes(label=sprintf("%.2f", 
                                                                                round(estimate, digits = 2))), nudge_y = .35, size=2.5)+
  theme_ggdist()+
  labs(x="",
       title="")+
  theme(legend.position = "none",
        axis.title.y = element_blank(),
        plot.title=element_text(color="black", face="bold", hjust=.5),
        axis.text.x=element_text(color="black"),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank())+
  annotate("segment", x = -.1, xend = -.07, y = .7, yend = .7, color = "black", linetype = "solid")+
  annotate("segment", x = -.1, xend = -.1, y = .7, yend = 1.3, color = "black", linetype = "solid")+
  annotate("segment", x = -.1, xend = -.07, y = 1.3, yend = 1.3, color = "black", linetype = "solid")+
  annotate("label", x = -.11, y = 1, label = "p<0.1", size = 2, color = "black")

plot_full+plot_UK+plot_ES+diff_full_plot+diffUK_plot+diffES_plot+
  plot_layout(ncol =3, nrow = 2, heights =c(3, 1))+
  plot_annotation(caption = 'Confidence intervals at 99%, 95% & 90% (N:6194)',
                  title = 'Study 1: Results from double-list experiment',
                  subtitle = 'Coefficients: Average treatment effect (ATE) of sensitive item') & 
  theme(plot.title = element_text(face="bold"))
ggsave("figures/Figure4.png", 
       height=6, width=10, dpi=600)
ggsave("figures/tifs/Figure4.tif", 
       height=6, width=10, dpi=600)


TableA3 <- list()
TableA3 [['All-Peers']] <-  lm_robust(list_answer ~ treat + round + list, data = df_het, clusters = IDvar)
TableA3 [['All-Self']] <-  lm_robust(list_answer_self  ~ treat + round + list, data = df_het, clusters = IDvar)
TableA3 [['UK-Peers']] <-  lm_robust(list_answer ~ treat + round + list, data = df_hetUK, clusters = IDvar)
TableA3 [['UK-Self']] <-  lm_robust(list_answer_self  ~ treat + round + list, data = df_hetUK, clusters = IDvar)
TableA3 [['Spain-Peers']] <-  lm_robust(list_answer ~ treat + round + list, data = df_hetES, clusters = IDvar)
TableA3 [['Spain-Self']] <-  lm_robust(list_answer_self  ~ treat + round + list, data = df_hetES, clusters = IDvar)
msummary(TableA3, star=TRUE, output = "tables/TableA3.tex")



df_het <- df_het %>%
  filter(gender != 2) %>%
  mutate(gender_bin = gender)
df_hetUK <- df_hetUK %>%
  filter(gender != 2) %>%
  mutate(gender_bin = gender)
df_hetES <- df_hetES %>%
  filter(gender != 2) %>%
  mutate(gender_bin = gender)
TableA4 <- list()
TableA4[['All-Peers']]  <- lm_robust(list_answer ~ treat*gender_bin + round + list, data = df_het, clusters = IDvar)
TableA4[['All-Self']]   <- lm_robust(list_answer_self ~ treat*gender_bin + round + list, data = df_het, clusters = IDvar)
TableA4[['UK-Peers']]   <- lm_robust(list_answer ~ treat*gender_bin + round + list, data = df_hetUK, clusters = IDvar)
TableA4[['UK-Self']]    <- lm_robust(list_answer_self ~ treat*gender_bin + round + list, data = df_hetUK, clusters = IDvar)
TableA4[['Spain-Peers']]<- lm_robust(list_answer ~ treat*gender_bin + round + list, data = df_hetES, clusters = IDvar)
TableA4[['Spain-Self']] <- lm_robust(list_answer_self ~ treat*gender_bin + round + list, data = df_hetES, clusters = IDvar)
msummary(TableA4, star = TRUE, output = "tables/TableA4.tex")


TableA5 <- list()
TableA5 [['All-Peers']] <-  lm_robust(list_answer ~ treat*factor(listround), data = df_het, clusters = IDvar)
TableA5 [['All-Self']] <-  lm_robust(list_answer_self  ~ treat*factor(listround), data = df_het, clusters = IDvar)
TableA5 [['UK-Peers']] <-  lm_robust(list_answer ~ treat*factor(listround), data = df_hetUK, clusters = IDvar)
TableA5 [['UK-Self']] <-  lm_robust(list_answer_self  ~ treat*factor(listround), data = df_hetUK, clusters = IDvar)
TableA5 [['Spain-Peers']] <-  lm_robust(list_answer ~ treat*factor(listround), data = df_hetES, clusters = IDvar)
TableA5 [['Spain-Self']] <-  lm_robust(list_answer_self  ~ treat*factor(listround), data = df_hetES, clusters = IDvar)
msummary(TableA5, star=TRUE, output = "tables/TableA5.tex")



